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Abstract 



In this paper, we extend our Smooth Particle Hydrodynamics (SPH) im- 
pact code to include the effect of porosity at a sub-resolution scale by adapting 
the so-called P — alpha model. Many small bodies in the different populations 
of asteroids and comets are believed to contain a high degree of porosity and 
the determination of both their coUisional evolution and the outcome of their 
disruption requires that the effect of porosity is taken into account in the com- 
putation of those processes. Here, we present our model and show how porosity 
interfaces with the clastic-pcrfcctly plastic material description and the brittle 
fracture model generally used to simulate the fragmentation of non-porous rocky 
bodies. We investigate various compaction models and discuss their suitability 
to simulate the compaction of (highly) porous material. Then, we perform sim- 
ple test cases where we compare results of the simulations to the theoretical 
solutions. We also present a Deep Impact-like simulation to show the effect 
of porosity on the outcome of an impact. Detailed validation tests will be pre- 
sented in a next paper by comparison with high- velocity laboratory experiments 
on porous materials (Jutzi et al., in preparation). Once validated at small scales, 
our new impact code can then be used at larger scales to study impacts and 
collisions involving brittle solids including porosity, such as the parent bodies 
of C-type asteroid families or cometary materials, both in the strength- and in 
the gravity-dominated regime. 
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1 Introduction 



The collisional process plays a key role in all stages of a planetary system history, 
from, the phase of planetary formation through collisional accretion to the late 
phases where the populations of small bodies evolve collisionally in a disruptive 
way. So far, these problems have been addressed using models appropriate to 
solid bodies represented as brittle rocky materials in which porosity was ne- 
glected or modeled at a macroscopic (i.e. resolved) level. Numerical codes, 
called hydrocodes (see, e.g., Benz and Asphaug, 1994) have been developed to 
compute the impact induced fragmentation of such solid bodies by solving the 
elastic- plastic conservation equations associated with a model of brittle failure to 
account for the fracture of the solid material. This already allowed to improve 
our understanding of the impact response of small bodies such as basalt-like 
material and asteroids belonging to the taxonomic class S, supposed to be com- 
posed of material with negligible porosity. However, several evidence point to 
the presence of a high degree of porosity in some small body classes. Asteroids 
belonging to the C taxonomic class are now believed to be highly porous, as 
indicated by the low bulk density (» 1.3 g/cm^) estimated for some of them, 
such as the asteroid 253 Mathilde encountered by the NEAR Shoemaker space- 
craft (Yeomans et al., 1997), and as inferred from meteorite analysis (Britt et 
al., 2006). Small body populations evolving at larger distances from the Sun, 
i.e. the Jupiter-Family Comets, the Kuiper Belt Objects, and the other classes 
of comets contain also a high level of porosity, as indicated by the estimated 
low bulk densities (below 1 g/cm^, e.g., Rickman, 1998) and the analysis of 
interplanetary dust particles collected on Earth. 

In parallel, the dissipative properties of porous media are invoked more and 
more as playing an important role in the formation of early planetesimals (e.g., 
Wurm et al., 2005). Hence, porosity emerges slowly as playing a major role 
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from the time of the formation of the planets to the coUisional evolution of the 
present day Solar System. 

Despite the growing focus on porosity, our ability to model its effect on the 
outcome of impacts and collisions remains limited. To address the question of 
modeling porosity in this context it is necessary to first define its scale in com- 
parison with the other relevant dimensions involved in the problem such as the 
size of the projectile and/or the crater, etc. In particular we define microscopic 
porosity as a type of porosity characterized by pores sufficiently small that one 
can reasonably assume that they are distributed uniformly and isotropically over 
these relevant scales. Macroscopic porosity on the other hand would be char- 
acterized by pores with sizes such that the medium can no longer be assumed 
to have homogeneous and isotropic characteristics over the scales of interest. 
In this case, pores have to be modeled explicitly and the current hydrocodes 
developed for the modeling of non-porous brittle solids can still be used. The 
presence of these large macroscopic voids will only affect the transfer efficiency 
and the geometry of the shock wave resulting from the impact which can be 
computed using this existing software. On the other hand, a body containing 
microporosity (i.e. porosity on a scale much smaller than can be numerically re- 
solved) can be crushable: cratering on a microporous asteroid might be an event 
involving compaction rather than ejection (Housen et al., 1999). In an impact 
in microporous material, a part of kinetic energy is dissipated by compaction 
which leads to less ejection and lower velocities of the ejected material. These 
effects cannot be reproduced using hydrocodes developed for the modeling of 
non-porous solids. Therefore, a model is needed which takes pore compaction 
into account. 

Sirono (2004) proposed a model which can be used to study low velocity 
collisions of porous aggregates. Using this model, the author showed that the 
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energy dissipation by compaction can lead to the sticking of dust aggregates. 
However, this model is not appropriate for high velocity impacts because of the 
simplicity of the equation of state (by construction, the pressure only depends 
on density). 

Recently, Wuennemann et al. (2006) proposed a so-called e — alpha model 
suitable to model porous material and for the use in hydrocodes. In this paper, 
we propose an alternative approach based on the so-called P — alpha model 
(Herrmann, 1969) which we found to be more appropriate to use with our nu- 
merical method (see Sec. 12. 2p . Our 3D Smooth Particle Hydrodynamics (SPH) 
code, originally developed by Benz and Asphaug (1994), includes a model of 
brittle failure of non-porous material and successfully reproduced impact exper- 
iments on non-porous basalt targets. Moreover, associated to the N-body code 
pkdgrav to account for the effect of gravity (see e.g., Richardson et al., 2000), 
it reproduced successfully for the first time the formation of S-type asteroid 
families resulting from the catastrophic disruption of non-porous parent bodies 
(see, e.g., Michel et al., 2001; Michel et al. 2003). A first implementation of a 
porosity model has been made in this code, by Benz and Jutzi (2006) using a 
simplified version of the P — alpha model (the density p was used instead of the 
pressure P) but was found to be inappropriate for materials with a high degree 
of porosity (see Sec. 12. 3p . Here, we improve on this first work by implementing 
the full P — alpha model which offers several advantages and has a more general 
application. 

In the following, we begin by describing our porosity model and comparing 
it to alternative models. We discuss some (counter intuitive) effects which occur 
when dealing with highly porous material. In the case of a high porosity, the 
energy dissipation during compaction leads to a large thermal pressure which can 
even become large enough to cause a decrease of the matrix density (Sec. 12. 3p . 
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We show that our model is suitable to model these effects. In Sec. [31 we recall 
the equations used in our modeling of brittle solids and then we show how 
porosity is related to these equations (Sec. S]). In Sec. [5l simple test cases are 
presented: first, a simulation of a ID shock wave in porous media is compared 
to the analytical solution; in the second test, we simulate the compression of 
porous pumice and measure the so-called crush-curve. As a first application, 
we perform a simulation of a Deep Impact-like impact and compare simulations 
with and without porosity model. We also compare different compaction models 
using a porous basalt target. Conclusions are then exposed in Sec. [51 A detailed 
comparison of numerical simulations with high-velocity impact experiments on 
porous material is presented in a forthcoming paper (Jutzi et al., in preparation). 



2 Porosity model 

While porosity at large scales can be modeled explicitly by introducing macro- 
scopic voids, porosity on a scale much smaller than the numerical resolution has 
to be modeled through a different approach. Our porosity model is based on the 
P — alpha model originally proposed by Herrmann (1969) and later modified 
by Carroll and Holt (1972). The model provides a description of microscopic 
porosity with pore-sizes beneath the spatial resolution of our numerical scheme 
(sub-resolution porosity) and which is homogeneous and isotropic. 

2.1 P — alpha model 

The basic idea underlying the P — alpha model consists in separating the volume 
change in a porous material into two parts: the pore collapse on one hand and 
the compression of the material composing the matrix on the other hand. This 
separation can be achieved by introducing the so-called distention parameter a 
defined as 
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where p is the bulk density of the porous material and ps is the density of 
the corresponding solid (matrix) material. Distention is related to porosity as 
1 — 1/a. According to its definition, the distention varies in the range ao > a > 
1, where is the initial distention. 

Using the distention parameter a, the equation of state (EOS) can be written 
in the general form: 

P^P{p,E,a) (2) 

According to CaroU and Holt (1972), the EOS of a porous material can explicitly 
be written as: 

P = -PsiPs,E,) = -P,{ap,E) (3) 
a a 

where Ps{ps,Es) represents the EOS of the solid phase of the material (the 
matrix). A crucial assumption in this model is that the pressure depends on 
the density of the matrix material. The pore space is modeled as empty voids 
and the internal energy E is assumed to be the same in the porous and the 
solid material {E — Eg), which implies that the surface energy of the pores is 
neglected. The factor 1/a in Eq. ^ was introduced by Carroll and Holt (1972) 
to take into account that the volume average of the stress in the matrix is given 

by 

Ps - aP (4) 

where P is the applied pressure. 

In the P — alpha model, the distention is solely a function of the pressure P: 



a = a{P) 



(5) 



where P = P{p,E,a). The relation between distention and pressure is often 
divided in an elastic regime {P < P^) and a plastic regime (P > Pe), where 
Pe is the pressure at which the transition between the two regimes occurs. In 
the elastic regime, the change of a with pressure is reversible. According to 
Herrmann (1969), the relation between distention and pressure can be defined 
as 



da 1 
dP Kq n(a) 



elastic - — [1 ~ (ut„:\2 )] (^) 



where Kq — CqPo and 



h{a)^l + ia-l) "° (7) 

co(ae - 1) 

where ae — a{P — Pe). This equation follows from the assumption that the 
elastic wave velocity changes linearly from the initial value Cg to the bulk sound 
speed Co in the solid state. Using Eq. ^ leads to a small change of the distention 
(from Q!o to Q!e) in the elastic regime and an elastic wave velocity which is smaller 
than in the case of a constant a. Since the distention changes only very little, 
it is often assumed that da / dP = in the elastic regime and therefore Ce = cq . 

In the plastic regime, the following quadratic form is often used to define 
the fmiction a = a{P): 

where Pe and Ps are constant. 

This is obviously a very simple model, but it is appropriate enough for 
many applications. A more realistic relation can be obtained experimentally 
by means of a one dimensional static compression of a sample, during which 
the actual distention am is measured as a function of the applied pressure Pm- 
The resulting crush-curve am{Pm) then provides the required relation between 
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distention and pressure for the material. In Sec. [51 we present test simulations 
where we use such a measured relation to define a(P). 



2.2 Alternative models 

According to the P~alpha model, distention is a function of pressure a — a{P). 
However, distention can also be defined as a function of other state variables. 
Wuenncmann et al. (2006) for example use the volumetric strain (e — alpha 
model). Another possibility is to define distention directly as a function of 
density. We will refer to that approach as p — alpha model. Note that in each 
model, the pressure is computed according to Eq. 

2.2.1 e — alpha versus p — alpha model 

According to Wuennemann et al. (2006), the volumetric strain can be expressed 
as 



where Vq is the initial volume and V the actual volume. Assuming that the vol- 
ume of the matrix is kept constant {Vs ~ Vso) volumetric strain and distention 
can be related as ey = ln{a/ao), which leads to the compaction function 



Another (very similar) way to define the evolution of the distention follows 
from its definition: a = ps/p. If we again assume a constant matrix volume and 
therefore a constant matrix density (ps — Pso), we get the compaction function 





a — aoe'^ 



(10) 




(11) 
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This relation also follows from Eqs. (|10p by replacing by 



Cv = ln{—) = ln{ — ) 



(12) 



where pg — Psq/chq is the initial density. 

Both equations, (I10|) and (fTT|) . describe the fasted possible way to decrease 
distention as a function of volume (density) change under the assumption that 
the matrix density is constant. However, as Wuennemann et al. (2006) point 
out, during the compression of a porous material, not only pore space is com- 
pacted, but also the matrix is slightly compressed leading to an increase of the 
pressure. Wuennemann et al. (2006) take this into account by introducing a 
parameter k to control the compaction rate: 



where eg is the critical strain where the compaction starts. Using Eq. and 
defining = ln{pQ/pe), we can write this equation in terms of density 



Wuennemann et al. (2006) also use a power-law compaction regime at a cer- 
tain volumetric strain. In this way, they are able to reproduce experimental 
compaction data obtained by a uniaxial compression of a porous material. In a 
similar way, one could also modify Eq. (jlip to obtain similar results. 

In both models, the time evolution of the distention parameter has a very 
simple form. In the e — alpha model, it is given by 



a = a^e' 



(13) 



(14) 



P 




da . 



(15) 
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and for the p — alpha model we get 



The comparison above shows that the e — alpha model and the p — alpha model 
are very similar and the only difference is the parameter which is chosen to 
measure the actual volume or density, respectively. 

2.2.2 p — alpha versus P — alpha model 

In the P — alpha model, the distention depends on the density via the pressure 
given by Eq. ([3]) but contrary to the two models described above, it is also a 
function of the internal energy. However, for small initial porosities, the energy 
contribution to the pressure remains small as long as a > 1, and can even 
be neglected in most cases. Therefore, the pressure can be approximated by 
P ~ P(p,a) and consequently, we can transform the function a{P[p,a\) in 
a — a{p). In this way, we can define a function a — a{p) which, instead of 
mimicking the behavior of the e — alpha model as in Eq. (|14p. approximately 
corresponds to the relation a ~ a{P), but neglects the thermal contribution to 
the pressure. Test simulations aimed at comparing a — a{P) and a — a{p) 
show that for low porosities this assumption is valid and the difference between 
the two models is rather small (see Sec. [5]). However, for high porosities this 
assumption is not valid and we observed some problems using this form of the 
p — alpha model. 

2.3 Problems with high porosities 

As we described above, the energy contribution to the pressure is very small in 
the porous regime (a > 1). However, this is only true for small initial porosities. 
In highly porous material, the thermal pressure can become large enough that 
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the density in the compressed state {a — 1) is below the initial density of the 
matrix (Zel'dovich and Raizer, 1967), while in a less porous material the density 
at a = 1 would be at least as high as the initial density of the matrix. Then, 
when the compressed state is reached, the density does not increase further with 
increasing pressure but rather decreases, and the volume increases accordingly. 
This behavior can lead to an anomalous but nonetheless actual Hugoniot curve 
(see Fig.[Tl top). It is important to stress that this anomalous behavior of highly 
porous material is not only a theoretical concept but can actually be observed 
in experiments. 

In the fully compressed state, by definition the material is only composed of 
the matrix. Since in this state the density of highly porous material is smaller 
than the initial value of the matrix density, this implies that the matrix den- 
sity in the fully compacted state has decreased from its initial value. Such a 
decrease must occur even before this fully compressed state has been reached 
(see Fig.[TJ bottom). As an important consequence, in highly porous materials 
the distention can be decreased faster (as a function of volume change) than in 
the case of a constant matrix density (Eqs. pU]) and pTjl ). This behavior oc- 
curs independently of the actual functional form of the distention. It is actually 
caused by the huge difference between the density of the porous material and 
the matrix density. 

It is important to note that even in highly porous material, the change of the 
matrix density during the compaction is very small compared to the change of 
the bulk density (see Fig. [1]). Nevertheless, it has a great effect since according 
to Eq. ([3]), the matrix density is used to compute the pressure. 

To correctly simulate the compaction of highly porous materials, therefore, 
the compaction functions must allow the density of the matrix material to de- 
crease. For the p — alpha or e — alpha models this implies that, above some 



13 



threshold pressure, the distension must decrease faster (as a function of volume 
change) than defined by Eq. (fTU]) or pT|) . which follow from the assumption of 
a constant matrix density. This could be achieved by a modification of Eq. (jl3p 
or (jl4p . On the other hand, defining distension as a function of pressure offers 
a simple prescription for compaction that allows for matrix expansion. The 
P — alpha model, without modification, can therefore simulate the compaction 
of highly porous materials given suitable model parameters. 

For illustration purposes, we show the pressure-distention relation in a unidi- 
mensional Shockwave under three different assumptions of compaction behavior, 
and for two initial porosities (Fig. In the first case we use the P — alpha 
model which assumes a quadratic function (Eq.[H]) for a = a{P) with Pe=8xl0® 
dyn/cm^ and ^5=7x10^ dyn/ cm^, and that ag = in the elastic regime. The 
values of Pe and Ps approximately correspond to the ones used by Herrmann 
(1969) to study a (low porosity) compaction wave in porous aluminium. For 
illustrative purposes, we assume that P^ and Ps do not depend on the initial 
porosity. In the second case we assume a constant matrix density during com- 
paction by using Eq. PH)) . with k = 1. The density is chosen so that the 
compaction starts at P = P^. As a third case, we also show the results obtained 
using a function a — a{p) which follows from the assumption that the thermal 
contribution to the pressure can be neglected and that the transformation of 
P ~ P(/0, a) holds, as described in the previous section. 

Since we only want to study the compaction behavior which takes place 
at moderate pressures {P < Ps) we use a simplified version of the Tillotson 
equation (Tillotson 1962, Melosh 1989) for this illustration. This simplified 
EOS provides a reasonable approximation of the full Tillotson equation in the 
considered regime. Furthermore, it allows us to examine the sensitivity of the 
pressure-distention relations on the equation of state (EOS) parameters. The 
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following equation is used: 



P = cpE + A/i 



(17) 



where c = a + b, fi = rj — 1 and rj = p/po and A, a and b are the usual Tillotson 
parameters. The assumptions made to obtain this equation are discussed in the 
appendix. 

The pressure-distention relation is finally computed using Eq. the sim- 
plified EOS (Eq. [TT)) and the Hugoniot equation 



where Vq — 1/po and V = 1/ p and Pq = and Eq ^ are used as initial values. 
We further use c = 2 and A — 7.5x10^^ dyn/cm^ (these values approximately 
correspond to the usual aluminium parameters) for the illustration (Fig. [2]). 

We have to point out the the following considerations are only valid in a 
moderate pressure regime where the assumptions leading to the simplified EOS 
are reasonable. 

We show two cases with a different initial distention, ao = 1.275 and = 
3.0. In the low porosity case, full compaction is reached with all three models at 
similar pressures. Of course, the curve obtained by using the P — alpha model 
corresponds to the quadratic relation used to define ol{P). The curve denoted 
by p — alpha shows the results obtained using the approximation P{p^ a, E) ^ 
P{p,a) followed by the transformation a{P) — a{P[p^a\) a — a{p), which 
assumes that the thermal contribution to the pressure can be neglected. As 
expected for low porosities, it has a similar shape as the P — alpha curve because 
in this case the thermal component of the pressure is small. The difference in 
pressure between the two curves is equivalent to the thermal contribution to the 
pressure at a given level of compaction. Note also that the P — alpha curve is 



E-Eo = {P + PoKVo-V)/2 



(18) 
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less steep than the constant matrix-density curve (e — alpha model curve with 
K — 1) because in this case the P — alpha relationship implies some compression 
of the matrix material during compaction. 

In the high porosity case, on the other hand, full compaction is achieved in 
the considered regime only with the P — alpha model. The curves of the other 
two models do not reach a — 1 until a much higher pressure is reached. Even 
the assumption of constant matrix density does not lead to full compaction. As 
described above, the difhculty in reaching complete compaction is caused by 
the fact that the matrix density can actually decrease with increasing pressure 
(due to the large thermal component) for highly porous materials. In the case 
of the constant matrix density, it can be shown that the minimal value that the 
distention can reach is given by 

(Xniin = lim a{e[P]) = uq — —^ ~ "0—7-;: (19) 

P^oo pq c + 2 c + 2 

We have to point out that this result follows using the simplified EOS (Eq. [T7)) 
which is of course not valid for infinite pressures. Using the full Tillotson equa- 
tion of state, the contribution of the term js^B and the reduction of the effective 
value of c at very high pressures (much higher than 10^*^ dyn/cm2) can still lead 
to full compaction. 

As Eq. shows, the parameter c determines whether or not full com- 
paction is possible at moderate pressures. For illustration, we once again show 
the curves which follow from the constant matrix-density assumption but now 
for different values of c (again for two cases with uq = 1.275 and uq = 3.0). 
For small porosities, changing c from 0.5 to 2 only slightly changes the slope 
of the curves (see Fig. [31 top). On the other hand, in the high porosity case, 
the parameter c has a great influence on the crushing behavior (see Fig. [31 bot- 
tom). For c = 0.5, full compaction is reached at P ~ 5x10^ dyn/cm^. However, 
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for c = 1 and c = 2 we get an amin of 1.0 and 1.5, respectively. The corre- 
sponding curves asymptotically approach these values. Again, using the full 
Tillotson equation, the behavior would be different for very high pressures and 
full compaction would be achieved even for c — 2. Nevertheless, we think that 
the estimation of amin (Eq. [TH]) can be used at least as a first guess of value of 
a critical distention 

c + 2 , , 

(XOcrit ^ (20) 

which follows from amm = 1- For an initial distention higher than this value 
(ao > Olocrit)^ auomalous effects can occur and models which do not allow matrix 
expansion fail to reach full compaction at moderate pressures. Obviously, this 
procedure to determine aocrit only works for EOS where a parameter c can 
be identified. In any case, the results indicate that using models where the 
distention is decoupled from pressure (e — alpha and p — alpha model), the 
compaction behavior can strongly depend on the EOS parameters which relate 
energy and pressure. On the other hand, using a model where the distention is 
a direct function of the pressure (P — alpha model), the compaction behavior 
is not sensitive to these parameters. 

In Sec. 15.3.11 results of impact simulations using the P — alpha, p — alpha 
and e — alpha model are compared. 

2.4 Our actual model 

Although the use of the p — alpha model (as an analogue to the e — alpha model) 
has some advantages, especially the simple form of the time evolution a — ^p, 
we found that for our numerical scheme, the most appropriate variable to define 
a functional form of the distention is pressure: a — a{P). The following reasons 
support this choice: 

• The P — alpha model can be used without modification to simulate the 
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compaction of highly porous material. 

• The relation between distention and pressure, which is used as an input 
in our model, can directly be obtained from the experimental crush-curve 
for the material considered. 

• We found that using our numerical scheme, no iteration is needed to im- 
plicitly solve for pressure and distention. 

For the actual form of the relation between distention and pressure we either 
use a quadratic relation ([5]) or, if available, we directly use the experimentally 
measured crush-curve to define the function a{P). By definition, the pressure 
distention relation {a{P)) we actually use in our code does not depend on the 
strain rate. Using a crush-curve which was measured under quasi static condi- 
tions therefore assumes that the same curve would be obtained at high strain 
rates. Experiments (Nakamura et al., in preparation) show that there actually 
is a (small) strain rate dependence of the crushcurve. However, we found (Jutzi 
et al., in preparation) that the results of impact simulations (i.e., the fragment 
mass distribution) do not strongly depend on the exact shape of the P — alpha 
relation. Therefore, we think that it is not problematic to use a low strain rate 
crush-curve as input in our code to simulate high strain rate events. 

The following function allows a good fit to a wide range of experimental 
crush-curves: 



a{P) 



(ae - at) lplZp})li + {at - 1) [plZp])!. + 1 ii P^ < P < Pt 

{at - ^) [p:Zp!)L +1 iiPt<P<Ps (21) 

1 if R < P 



where Ps,Pe and ae have the same meaning as in Eq. ([H]), and P^ < Pt < Ps 
and 1 < at < ao are parameters indicating a transition pressure and distention. 
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respectively. The function (j2ip and its first derivative are smooth by definition, 
which aUows the existence of two regimes of a{P), each with a individual slope 
(ni and 71.2). 

In the elastic region {P < P^) we either use Eq. © to define [da/dP]eiastic 
or, to simplify matters, we assume that the distention is constant, i.e. a,, — ag. 
For all simulations presented in this paper (except the one in Sec. 15. ip . we use 
[da/dP]eiastic = in the elastic regime. 

The derivative da/dP which is used to compute the time evolution of the 
distention is now given by 

[da/dP]elasUc liP <Pe 

da/dP = { (22) 
[da/dP]piasUc otherwise 

where [da / dP]piastic follows from Eq. (|2ip. We further assume that unloading 
(from a partially compacted) state is elastic. Consequently, we define 

da/dP iidP>0 
da/dP = { (23) 

[da/dP]eiasUc otherwise. 

As discussed above, [da/dP]eiastic is either computed using Eq. ^ or assumed 
to be zero. 

The time evolution of the distention parameter can be written as 



da dP 

a = (24) 

dP dt ^ ^ 



Using Eq. ([3]) we finally get 



e{§^)+-p{W:) da 



«w = \:l • ^ (25) 



da 
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The equations ^ and (^11 - define the constitutive equation which de- 
scribes the compaction behavior of a porous material. In the original work of 
Herrmann (1969) and Carroll and Holt (1972), the P — alpha model was in- 
tended to be a first order theory in which shear strength effects are considered 
secondary and consequently, the stress tensor was assumed to be diagonal. In 
this work, we use the full stress tensor and therefore, we extend the original 
model with a relation between the distention and the deviatoric stress tensor 
(Sec. O and 1121). 



3 Model equations and implementation 

Our numerical technique is based on the Lagrangian Smooth Particle Hydro- 
dynamic (SPH) method. Since the basic method has already been described 
in many papers (see for examples reviews by Benz, 1990; Monaghan, 1992) we 
refer the interested reader to these earlier papers. 

The standard gas dynamics SPH approach was extended (see for example 
Libersky and Petschek, 1991) to include an elastic-perfectly plastic material 
description and a fracture model based on the one of Grady and Kipp (1980) in 
order to model the behavior of brittle solids (Benz and Asphaug, 1994). As our 
porosity model interfaces with this material description, we begin with a short 
review of this previous approach. Note that the following equations describe 
non-porous material. Porosity is introduced in Sec. 21 

3.1 Elastic perfectly plastic strength model 

The equations to be solved are the well-known conservation equations of elasto- 
dynamics; they can be found in most standard textbooks. The mass conserva- 
tion can be written as: 



dp" dv'^^ 
dt ^ dx 



PlwT^O (26) 
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where dj dt is the Lagrangian time derivative, p the density, v the velocity and 
X the position. The conservation of momentum has the following form: 



dv'^ _ 1 da'^^ 

~dt~~p~dx>^ ^ ' 

where ct**^ is the stress tensor given by 

^ gK\ _ p^K\ ^28) 

where P is the hydrostatic pressure, 5'^^ is the Kroneker symbol and S*^^ is the 
(traceless) deviatoric stress tensor. Finally, the conservation of energy is given 
by the equation 

^ = _:^J_^« + l5-A..A (29) 
dt p ox'^ p 

where e is the strain rate tensor given by 



" 2 [dx^ ~^ dx^ 



(30) 



In order to specify the time evolution of the deviatoric stress tensor S'^^ we 
adopt Hooke's law and define the time evolution of the deviatoric stress tensor 
as: 

— — = e"^ - S'^^e"" + S'^^n^" + S^^n"" (31) 
dt V 3 / 

where fj, is the shear modulus, and Cl is the rotation rate tensor: 



""'^rs-sii^ (32) 



Finally, plasticity is treated using the von Mises yielding criterion. 

In order to solve this set of equations, an equation of state has to be specified 
which relates density, energy and pressure: 
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P = P{p,E) 



(33) 



For the simulations presented in this paper we use the so-called Tillotson equa- 
tion of state (e.g., Tillotson 1962, Melosh 1989). 

3.2 Fracture 

Brittle materials cannot be modeled using elasticity and plasticity alone because 
these materials fail under tension or shear stress. To take this behavior into ac- 
count, wc use a fracture model introduced by Grady and Kipp (1980) and based 
on the presence of incipient flaws in the material and on crack propagation under 
increasing strain. This model has been introduced using an explicit distribution 
of incipient flaws by Benz and Asphaug (1994, 1995) in their SPH code. In this 
model it is assumed that the number density of active flaws at strain e is given 
by a Weibull distribution (WeibuU, 1939) 

n(e) = fee™ (34) 

where k and m are the material dependent Weibull parameters. When the local 
tensile strain has reached the activation threshold of a flaw, a crack is allowed to 
grow at a constant velocity Cg which is some fraction of the local sound speed. 

The half length of a growing crack is therefore 

a = Cg{t-t') (35) 

where t' is the crack activation time. 

Crack growth leads to a release of local stresses. To model this behavior, 
we follow Benz and Asphaug (1994, 1995) and introduce a state variable D (for 
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damage) which expresses the reduction in strength under tensile loading: 



ao - fT(l - D) (36) 

where a is the elastic stress in the absence of damage and (Td is the damage- 
relieved stress. The state variable D is defined locally as the fractional volume 
that is relieved of stress by local growing cracks 

= 3_ (37) 

where V = A/iirR^ is the volume in which a crack of half - length is growing. 
Using Eqs. ([55)1 and ([57]) we get the following equation for the damage growth 

Damage accumulates at a rate given by Eq. (|38p when the local tensile strain 
€i reached the activation threshold of a flaw. Note that Ci is obtained from the 
maximum tensile stress aj after a principal axis transformation 



(39) 



(l-A)-B 

where Di is the local value of the damage and E is the Young modulus. 

4 Interfacing porosity with the sohd model 

So far, we only described how porosity (i.e. the distention parameter a) is used 
to modify the pressure: 



a 
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In this section, we show how distention interfaces with the material model ex- 
posed in the previous section, which is used to describe the behavior of solids 
under strain increase. 



4.1 Distention and strength 

As we have discussed in Sec. [2l the pressure P is calculated using the matrix 
density ps instead of p. Consequently, the deviatoric stress tensor has to be 
modified as well. In order to compute the time evolution of S"^^ as a function 
of the matrix variables, we introduce the following factor: 



This factor relates the velocity divergence of the matrix and of the porous ma- 
terial. Using the continuity equation for the matrix 



Ps = -ps[^v\s 



(41) 



and for the porous material 



p[W] 



(42) 



we can write the factor / as 




Ps_P _ Ps_ 

ps P ap 



(43) 



Using Ps ~ ap + ap wc finally get 




(44) 
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The factor / is then used to compute the time evolution of S'^^ for the porous 
material: 



dt 



f 



dS"^ 
dt 



(45) 



The multiplication by the factor / is motivated by the fact that both, the 
velocity divergence 



Vv 



dvi dv2 



dv?. 



(46) 



dxi 0x2 dxs 

and the time derivative of the deviatoric stress tensor (Eq. I3ip are linear com- 
binations of the spatial derivative of the components of the velocity vector (the 
linearity of Eq. pTjl follows from Hooke's law). Since according to Eq. the 
velocity divergence of the matrix is given by 



(47) 



we obtain 



dS'' 



dt 



= f 



dt 



(48) 



In addition to the multiplication by /, the deviatoric stress tensor S'"'*' is mul- 
tiplied by a^^ as it is done with the hydrostatic pressure P. We finally write 
the time evolution of S'^^ in the following form: 



d 
di 



1 



-S 



kA 



1 dS""^ 1 cxda 



dt 



dt 



(49) 



where dS'^^ /dt is modified according to Eq. (|45p . 

The computation of the factor / can fail for small p. This is the main reason 
why the p — alpha model was used in our first implementation (Benz and Jutzi, 
2006) as in this case, this factor can be computed using a simpler relation: 
/ = {p/a){da/dp). For several reasons (see Sec. 2.4), the P — alpha model is 
actually more appropriate. Therefore, we have worked out a functional form for 
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the factor / that does not lead to difficulties for small p. For this, we replace E 
in Eq. ([^ hy E — P/ ■ p and we rewrite Eq. (|^ as 



Defining 



a - 



da 
dP 



P 



da dp 
dP ' 'dt 



da^ P/p'{^)+-{^) 
dp 



da 
dP 



the derivative of a can be written as 



da 
dP 



(50) 



(51) 



da . 

a = —p 
dp 



(52) 



and we finally compute the correction factor in the following form: 



/ = 1 + ^ = 1 + ^^ 
ap dp a 



(53) 



The time evolution of the deviatoric stress tensor is then computed using Eq. (|49|) 
and (|53p . In this way, the deviatoric stress is a function of the distention and 
and also of the used P — alpha relation. On the other hand, we do not explic- 
itly relate the yield strength Y (used for the von Mises yielding criterion) and 
distention. 



4.2 Distention and damage 

Porosity does not only affect the stress behavior. It also has to be taken into 
account to compute the state variable damage. 

Compression of a porous material is accompanied, if significant enough, by 
the breaking of cell walls. Our model takes into account this crushing behavior 
relating distention with the state variable damage. Since both damage D and 
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distention a are defined as volume ratios (Eqs. ([57)1 and ([T]), respectively), we 
assume for simplicity a linear relation between D and a (other forms will be 
investigated in the future). The conditions Z) = at a = ao, and Z? = 1 when 
all pores have been crushed {a = 1), lead to the following expression: 



D = 1 



The time evolution of D^l^(ci) is then given by 



dt 



da dt 



(54) 



(55) 



and using Eq. ([51]) we obtain 



dt ^ ^3 



1 - 



a-1 



QfO - 1 



1 c?ck 



1 dt 



(56) 



A close examination of Eq. (j56p reveals a problem since for a = ao, the derivative 
D^/^ /dt becomes infinite. In order to avoid this problem we add the small 
quantity 5D to the linear relation l|54p and normalize it so that the conditions 
Z3(ao) = and D{1) ~ 1 are still satisfied: 



D 



1/3 



{D + 5Df''^ - (SD) 



Nl/3 



(l + 5£))l/3 

Using Eqs. ([55)) and we finally get 



((5L')i/3 



(57) 



dt 



1 _ + 
(qo-i) 



da 



3 (1 + ^L*)^/-'^ - (,5L>)i/3 (ao - 1) dt 



(58) 



where we set SD — 0.01. The actual relation between distention and damage is 
shown in Figure [H 

We now have two equations describing damage growth: the first treats dam- 
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age under tension while the second ([5^ is related to the compression of 
the (porous) material. Note that we do not explicitly model damage increase 
due to shear deformation. However, since the local scalar strain is computed 
from the maximum negative stress after principal axis transformation (Eq. I39p , 
shear fracturing is implicitly accounted for in our model. 

In order to get the total growth of damage, we build the sum of the two 
differential equations: 











-dD^/^' 


dt 


total 


dt 


+ 

tension 


dt 



(59) 



We now use this equation instead of Eq. (j38|) to compute damage D which varies 
between and 1. According to Eq. p6|) . damage leads to a reduction of both 
tensile and shear strength by a factor of (1 — D). 

According to our model, damage can grow only. We do not include any 
restoration of damage as Sirono (2004) proposed for low density grain aggre- 
gates. 

4.3 Material parameters 

All parameters used by our porosity model are material parameters which can 
in principle be measured experimentally. Some of these parameters, such as 
the crush-curve, can be measured more easily than others (such as WeibuU 
parameters, shear and tensile strengths). Unfortunately, such measurements 
are rarely done in practice and we plan to carry out some of them for several 
porous materials in future works. 

The lack of an experimentally determined reliable database of relevant ma- 
terial parameters is actually one of the most limiting factor in our model. In 
particular, the thorough testing of the model by comparison with experiments 
is rendered particularly difficult if all material properties have not been mea- 
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sured properly. Freely choosing the missing values so as to match an experiment 
is not a satisfactory approach for an ab initio method such as ours. Unfortu- 
nately, this is often the only alternative we have. However, recent experiments 
have been performed on some porous materials, and material properties have 
been measured. A detailed comparison between results of simulations using our 
model with these experiments will be the subject of a next publication. 

5 Tests 

In this section we present some simple test cases where we compare our numeri- 
cal model to expected theoretical solutions. We also present a Deep Impact-like 
simulation to show the effect of porosity on the outcome of an impact. Further, 
different compaction models are compared by performing an impact in (highly) 
porous basalt. 

The present tests aim at a first actual verification of the model by showing 
that it is consistent and correctly implemented in our code. A more detailed 
validation by comparison with actual impact experiments on porous material 
will be presented in a forthcoming paper. 

5.1 ID compaction wave 

As a first numerical test we consider a (plane) shock wave travelling in only 
one spatial dimension. We compare the simulation with the analytical results 
obtained by solving the corresponding Hugoniot equations. 

In this simulation, we use porous aluminium with an initial distention of 
ao — 1.275. We use a P — alpha relation with both an elastic and a plastic 
regime. In the elastic regime P < Pe, the distention is computed using Eq. ([6]). 
Using this equation instead of da/dP — leads to a smaller velocity of the elastic 
wave. We actually use Eq. ([6]) in this case because it provides an additional test 
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(properties of the elastic wave). However, for all other simulations presented in 
this paper we assume that — and therefore Cg = cq. 

For the plastic regime we use the quadratic relation ([H]). The porosity pa- 
rameters used in this simulation are given in Table [T] 

As for the equation of state, we use the (full) Tillotson equation with alu- 
minum parameters (Melosh, 1989). Since we are dealing with a purely one 
dimensional problem, strength is not included in this simulation and only the 
diagonal part of the stress tensor (pressure) is considered. 

The analytical solution of this problem can be obtained by solving the Hugo- 
niot equations together with Eq. ([3]) and the P — alpha relation defined above. 

To carry out the simulation we use a cylinder aligned along the z-axis with 
a radius r — 0.2 cm and a height h — 2 cm. Since we want to study a one 
dimensional case, the forces acting on the particles in the x-y plane are set to 
zero: 



The shock wave is then produced by moving the particles of the first layer (at 
the top of the cylinder) with a constant velocity Vz = -45.8 x 10"^ cm/s in 
the z-direction. The last layer (at the bottom of the cylinder) is fixed. We 
use 5.6 X 10^ particles in this simulation. However, in order to reduce boundary 
effects, only particles in a cylinder of 0.05 cm radius are used for the comparison 
with the theoretical solution. 

Figure [5] (top) shows the compaction wave travelling in the z-direction at 
time t = 3.5 fjs. As expected, there are two waves: first, there is a so-called 
elastic precursor with an amplitude equal to the elastic pressure P^. The elastic 
precursor is followed by a plastic compaction (shock) wave. The distention is 
only slightly changed due to the elastic wave and is decreased to a = 1 by the 
compaction wave (see Fig. [SJ bottom). 
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There is a good agreement between the theoretical solution (solid hne) and 
the velocities of the two waves (the elastic precursor and the shock wave) ob- 
tained by the simulation. Moreover the corresponding pressure amplitudes have 
the expected value. The oscillations behind the compaction wave are due to the 
moving layer of particles (wall). The smoothed shape of the wave (especially of 
the elastic precursor) is caused by the artificial viscosity and the smoothing due 
to the SPH interpolation. No effort was made to fine tune the artificial viscosity 
to achieve a less diffusive solution. 

5.2 Crush-curve simulation 

As we have discussed in Sec. [U an advantage of defining the distention as a 
function of pressure is that this relation can be directly determined experimen- 
tally by compressing a porous sample. Such an experiment (Nakamura et al., 
in preparation) was performed by A.M. Nakamura and K. Hiraoka at Kobe 
University (Japan). In this experiment, a cylindrical sample of porous pumice, 
confined within a steel cylinder, is compressed in one dimension. The applied 
force F and the corresponding displacement e of the penetrating piston (which 
defines the actual length / of the sample) are measured. From these quantities, 
the applied pressure P„i is given by the force F per unit area and the actual 
distention can be determined in the following way 

I , . 

OLm = OLQ— (60) 

to 

were is the initial length of the sample and ^ = ^ e is the actual length. 
In this way, we can obtain the relation OLmiPin)- Using the function (|2ip . we 
fit the curve OLmiPm) and get the required analytical relation a = a{P). The 
fitting parameters are given in Table [H 

As a further test, we use this relation in our porosity model and we simulate 
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the compaction experiment. The simulation is performed moving the first layer 
of particles of a cylinder A aligned along the z-axis with the velocity v^. As in 
the experiment, the walls and the bottom of the cylinder are fixed. In order 
to determine the actual (macroscopic) distention and the applied pressure, we 
define a test cylinder B within A with a radius r and an initial length of 
corresponding to the sample in the experiment. Since it is a quasi-static com- 
pression, we assume that the applied pressure corresponds to the pressure in 
the sample and we compute P as the average pressure in the cylinder B. As in 
the experiment, we define the actual distention using equation Eq. (j60|) . 

Figure [H] illustrates the setup of the simulation which has been performed 
in 3D. On this 2D slice, the dark particles correspond to the cylinder B which 
represents the sample. The simulation is shown at the initial state (top) and 
during the compression (bottom). In order to control the uniformity of the 
compaction, we compute the standard deviation of the distention in the sample 
(particles within B) at certain timesteps. During the whole simulation, the 
standard deviation is of the order of 1% which indicates that the compaction 
proceeds uniformly. 

In Figure [71 the crush-curve obtained by the simulation is compared to the 
experimental crush-curve. There is a very good agreement between the two 
curves. Of course, using the measured crush-curve as an input, one expects to 
obtain the same curve by the simulation. However, we have to point out that 
the pressure - distention relation used as input in the code is not necessarily 
the same as the pressure - distention relation measured in the simulation. In 
the a{P) relation used in the code, distention is defined as a = Ps/p- In the 
relation obtained by the simulation, distention is defined according to Eq. (|60p , 
which can be written as: 

otm = aay- = — (61) 

'0 PsO 
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This definition therefore assumes a constant matrix density {ps — Pso)- As we 
discussed in Sec. 12.31 the matrix density changes during the compaction. How- 
ever, the change of ps is only small (typically 1%) leading to an error of a„i 
which is in the same order. The fact that the two curves (simulated crush-curve 
and the crush-curve used as input) are in a very good agreement indicates that 
the change of the matrix denstiy is indeed very small. This result therefore sug- 
gests that the distention defined by am — Pso/p is a reasonable approximation 
of the true distention a — Ps/ p (even for highly porous material). However, it is 
important to note that even though the matrix density changes only very little, 
it can still have a great effect (see Sec. 12. 3p . 

In a quasi static compression, the energy increase {PdV work) is generally 
lower than in the case of shock compression, leading to a smaller thermal pres- 
sure. In the simulation above we assumed that the system is thermally insulated 
(i.e., the energy radiation is not taken into account). Therefore, the thermal 
pressure can still not be ignored. In Fig. [71 the crush-curve obtained by a 
simulation using a relation a = a{p) which follows from the a — a{P) curve 
neglecting the thermal pressure {p — alpha model) is shown. To reach a given 
distention, a higher pressure is needed using the p — alpha model instead of the 
P — alpha one. This difference corresponds to the thermal pressure. 

5.3 Impact simulations 

In this section, we show several impact simulations using different material types 
(ice and basalt) and compaction models (no porosity model; P~alpha, p~alpha 
and e — alpha model). Note that the following results were obtained using a pre- 
damaged (i.e. strengthless) material. Therefore, we are in the gravity regime 
where the final crater size and the total amount of ejected material depend on the 
gravity acceleration. Due to the high resolution and the resulting small timestep 
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(which is much smaller than the crater formation time), we only simulate the 
first phase of the crater formation. 

5.3.1 Deep-impact like impact 

As a first application of our model, we show the effect of porosity on the out- 
come of an impact event by comparing simulations using porous and non-porous 
targets. Since we want to model a realistic case, we simulate a Deep-Impact- 
like impact, inspired from the Deep Impact Space mission (A'Hearn and Combi, 
2007). 

We model the target in two different ways. In the first case, we use an initial 
distention of ao = 3.0 and we explicitly model porosity. To simplify matters we 
use a quadratic relation (Eq. [8]) for a — a(P) with = fx 10^ dyn/cm^ and 
Ps = 2x10^ dyn/cm^. These values are chosen rather arbitrarily, however, the 
resulting a{P) relation looks reasonable compared to measured crush-curves of 
other materials (see for expample Fig. [7|). In the second case, we use the same 
initial density (po — Pso/oio)^ but the target is modeled as a solid (without 
porosity model). As target material we use pre-damaged (strengthless) ice. 
Only a small part of the target (comet Tempel-1) is modeled (half sphere with a 
radius of 28 m) . In order to have a reasonable resolution in the region of interest 
(crater), we use a rather high number of particles for the target {Nt — 5.2x10®). 
The resulting mass per SPH particle is mp — 2.6 kg and the smoothing length is 
= 23 cm. The impactor is modeled as a 370 kg aluminium sphere impacting 
at an angle of 30 degrees (from horizontal) with a velocity of 10 km/s. We use 
140 SPH particles for the impactor in order to have the same mass per particle 
as in the target. The Tillotson equation of state (without simplifications) is 
used for these simulations. 

Figure [5] shows the outcome of the simulation after 50 ms in two dimensional 
slices of the three dimensional target. The dark particles in these plots have 
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vertical velocities greater than 2 m/s (which is about the escape velocity of 
Tempel-1). Obviously, there is much less material ejected in the case where we 
explicitly model porosity (top) than in the solid target simulation (bottom). In 
fact, the difference of ejected mass is about a factor ten. Figure [5] shows the 
amount of mass ejected with a velocity higher than a certain velocity in both 
cases. 

For comparison, the results obtained by using the p— alpha model (instead of 
P ~ alpha) are also plotted in Fig. [5) For this, we use a relation a — a{p) which 
follows from a = a{P) by neglecting the thermal component of the pressure 
(see Sec. 12.2.2)) . There is only a very small difference between the results of the 
two models. This is consistent with the estimation fEa. l^O]) of aocrit — 6 for c = 
0.4 (ice). Since ao = 3 < aocrit, full compaction is possible and since c is very 
small, we do not expect big differences. However, the situation looks different 
if we use basalt instead of ice as target material. 

5.3.2 Comparison between different compaction models 

For the typically used basalt parameters we get c ~ 2 and therefore aocrit — 2. 
Consequentely, full compaction (at moderate pressures) is not possible for ao = 
3 > aocrit using a compaction model which does not allow an expansion of the 
matrix. 

To illustrate the different compaction behaviors using different compaction 
models, we simulate an impact in a basalt target with ao = 3.0. For this 
simulation, we use the same initial conditions as above but we use a smaller 
target and the impact is now head on. We investigate three cases using the 
follwing compaction models: 

1. P ~ alpha with the same parameters as above. 

2. p — alpha, using the relation a — a{p) which follows from a = a(P) 
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neglecting the thermal pressure 

3. e — alpha, k — 1 (constant matrix density assumption) 

In Fig.[Tni the degree of compaction (i.e., the distention a) of the target material 
due to the impact is shown for simulations using the model 1 (top), model 2 
(middle) and model 3 (bottom). On these plots, the black color corresponds to 
a distention of a — ao and the white color to full compaction a = 1. The results 
are shown at a time (20 ms) where compaction is finished and the distention a 
does not change anymore. 

As it can be seen, only in the simulation using model 1 there is a small fully 
compacted zone below the crater. Fig. [TT] shows the corresponding distention 
profiles. The minimal distention reached by model 2 and 3 is in both cases 
about a — 1.8. This indicates that also the maximal density is about the same 
in both cases, because a — Ps/ p and ps changes only very little (model 2) or is 
even constant (model 3). The total volume which is compacted (derived from 
the actual distention) is 62.3/61.9/82.9 m^ for the models 1/2/3. This means 
that even though model 3 does not lead to full compaction, the total volume 
of compacted material is higher than using model 1. The reason is that the 
distention profile resulting from model 3 is much less steep than using model 1 

In Fig. [T^ the cumulative volume of ejected material as a function of the 
vertical velocity is shown. The horizontal lines in the plot represent the total 
compacted volume. Interestingly, using model 1, there is less ejection than 
using model 3 (even though model 3 leads to more compaction). This can be 
explained by the fact that a given pressure amplitude leads to a higher degree 
of compaction using model 1 than using model 3, and therefore, more energy is 
dissipated. Using model 2, we get the same amount of ejected volume at high 
velocities as with model 3, but there is more material ejected at low velocities. 



36 



In order to study the dependence of the compacted and ejected volume on 
the compaction parameters, we perform a simulation using model 1 with Pg = 
1x10^ dyn/cm^ instead of Pg = 2x10^ dyn/cm^. The other parameters are the 
same as above. Fig. ll3l shows the compacted volume and the cumulative volume 
of ejected material for the two simulations using model 1 (each with a different 
Ps). Clearly, using a very low value for Pg (1x10^ dyn/cm^) leads to much 
more compaction (95.2 m^) and less ejection than using Pg = 2x10^ dyn/cm^. 
There is even more compaction than from model 3 which assumes a constant 
matrix density. 

6 Conclusions 

In this paper, we have presented a new approach to model sub-resolution poros- 
ity in brittle solids that can be coupled to a 3D SPH hydrocode in order to sim- 
ulate impacts and collisions involving porous bodies. Such bodies are believed 
to be present in all the populations of small bodies in our Solar System. There- 
fore, understanding their impact response is crucial to determine the collisional 
evolution of those populations during all stages, and to assess the accretion ef- 
ficiency of small bodies during planetary formation. Moreover, this can help 
defining efficient mitigation strategies against the impact of a porous body on 
Earth. 

In practice, the implementation of our model does not consume excessive 
CPU time and is easily implemented in a parallel code (porosity is a local 
property) so that simulations involving multi-million particles can readily be 
performed. In fact, extensive testing has shown that high resolution is really 
needed to obtain converged solutions in the case of simulations involving frac- 
turing and/or porosity. 

We presented two test cases. In the first test, we compared the simulation 
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of a ID compaction wave in porous material with the theoretical solution given 
by the Hugoniot equations. The amplitude and velocity of the resulting waves 
(elastic precursor and shock) in the simulation agree with the theoretical solu- 
tion. The second test is a simulation of an experiment consisting in compressing 
a porous sample. We show that an advantage of the P — alpha model is that an 
experimentally measured crush-curve can be directly used to define the relation 
a{P). Therefore, if a crush-curve has been measured for a given material, we 
can then simulate an impact on this material using its very crush-curve and not 
an arbitrary one which may lead to different outcomes. 

As an application of our model, we presented the simulation of a Deep 
Impact-like impact. The main conclusion of this first application is that to 
model the behavior of a porous material during an impact, it is not enough 
to use the small bulk density expected for a porous material and compute its 
response from an impact by using a classical model of brittle failure of non- 
porous rock material, as one might be tempted to do. Actually, the cratering or 
disruption of porous material involves different processes than the ones involved 
in a non-porous brittle material (e.g. the crushing of pores), and this makes a 
huge difference in the outcome. 

We also investigated the compaction behavior in an impact in basalt using 
different compaction models. We showed that using the p — alpha or e — alpha 
models (without modifications) can lead to difficulties to reach full compaction 
(a = 1) in the case of highly porous material. Using these models, the com- 
paction behavior depends much on the EOS parameters which relate energy and 
pressure. On the other hand, the P— alpha model can be used without modifica- 
tions to simulate the compaction of highly porous material and the compaction 
behavior is not sensitive to these EOS parameters. 

Our next step will be to validate our code in more details and in a context 
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adapted to its future applications by confronting it with real high-velocity im- 
pact experiments on porous targets. Such a validation will mainly consist of 
reproducing the size distribution of the fragments (and their velocities when 
they are measured), using as inputs all measured material parameters. Indeed, 
it is important that matching the results does not rely on a fine tuning of all 
parameters but rather on the validity of the physical model. This will be the 
subject of a next paper. Comparisons with impact experiments in laboratory 
will have to be performed in the long run for a wide variety of materials (porous 
and non porous) to improve our understanding of the impact process as a func- 
tion of material properties and impact conditions. The formation and evolution 
of planetary systems involves bodies with a wide range of properties and col- 
liding with each other in different regimes of impact energies, leading either to 
their accretion and the formation of planets, or to their disruption as in the 
current stage of our Solar System. It is then important to have the possibility 
to simulate the collisional process in these different regimes and between bod- 
ies with different degrees of porosity. Collisional evolution models, which need 
constraints to characterize the outcomes of collisions, will certainly benefit from 
these investigations. 
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Appendix 



Simplified Tillotson EOS 

First, we assume that the energy involved does not exceed the energy of incipient 
vaporization {E < Eiy). In this state, the EOS has the following form: 



pE + Afi + Bn^ (62) 



E/ [E^T^^] + 1 



where fi = rj ~ I and 77 = p/ Po and a, 6, A and B are Tillotson parameters. We 
further assume that /i remains small and therefore /i^ << fi. This assumption 
is motivated by the fact that according to Eq. ([3]) /i is given by /i = 1 — PsI PsO 
and even though the matrix density is not constant (as discussed above), the 
variation remains small as long P < Ps (about 1% in Fig. [T]). In fact, it can be 
shown that if Ps << A then (ps — Pso)/psO << 1- According to Eq. the 
pressure is computed using the matrix density (p = 1 — ps/ pso), and therefore 
/i^ ~ is a reasonable assumption. 

Finally, we only want to consider cases where the energy E remains small 
compared to the parameter Eq. The energy at P = can be estimated by: 

2 pq p 2 pso pso 2 pso 
and therefore it is required that 

E^IPs^^«Eo. (63) 
2 PsO 

For Ps— 7x10^ dyn/cm^, pso—2.7 g/cm^ and ao=3.0 we get E ~ 2.6x10^ erg/g 
which is indeed small compared to a typical value of Eq (e.g. 5x10^° erg/g for 
aluminium) and so this condition is fullfilled. 

Using the above assumptions, we can rewrite Eq. (|62p in the very simple 
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form 



P = cpE + A^i 



(64) 



where c = a + 6 ~ constant. 
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Pe 


8e8 dyn/cm'^ 




7e9 dyn/cm^ 


Co 


5.35e5 cm/s 




4.11e5 cm/s 



Table 1: Parameters used in our porosity model for porous aluminium with an 
initial distention ao = 1.275. Definitions of parameters are given in the text. 
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ao 


4.64 


at 


1.90 


Pe 


1.00e7 dyn/cm^ 


Pt 


6.80e9 dyn/cm^ 


Ps 


2.13e9 dyn/cm^ 



Table 2: Parameters used to fit the crush-curve of pumice. 
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Figure 1: Pressure versus density (top) and pressure versus matrix- density 
(bottom) for a material with low ( ao = 1.275) and high (ao = 3.0) porosity. 
In the case of a high porosity, the density does never reach the initial matrix 
density (dashed horizontal line) and the matrix density itself decreases at a 
certain point, before a = 1 (dashed vertical line). 

Figure 2: Comparison of the pressure-distention relationship for different com- 
paction models. The P — alpha model assumes a quadratic function for a = 
a(P). For the p — alpha model we use a function a = a{p) which follows from 
a = a{P) by neglecting the thermal pressure. In the e— alpha model we assume a 
constant matrix density (k = 1). In the low porosity case (top), full compaction 
is reached with all three models. In the high porosity case (bottom), only the 
P — alpha model leads to full compaction. The minimal distention obtained 
using the constant matrix density assumption is about a=1.5. 

Figure 3: Pressure-distention relationship for different values of c ~ a + b 
(Tillotson parameters) using the constant matrix density assumption. For low 
porosities, the influence of c is rather small. In the high velocity case, the value 
of c determines whether or not full compaction is possible at moderate pressures. 

Figure 4: Relation between distention and damage. 

Figure 5: SPH Simulation of a ID compaction wave in a 3D cylinder composed 
of porous aluminium. There are two waves (top): an elastic precursor followed 
by the compaction (shock) wave. Bottom: the main decrease of the distention 
is due to the compaction wave. 

Figure 6: Simulation of compression of a porous sample shown when a = ao = 
4.64 (top) and during the compression (bottom) when a = 2.2. The sample is 
represented by the dark particles. 
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Figure 7: Crush-curve of porous pumice measured by A.M. Nakamura and K. 
Hiraoka at Kobe University and obtained by a simulation using the P — alpha 
and p — alpha model. 

Figure 8: Simulation of a Deep Impact-like impact in porous ice with ao = 3.0 
(top) and non porous ice with the same initial density (bottom). Dark particles 
have a vertical velocity greater than 2 m/s. 

Figure 9: Cumulated mass of ejecta as a function of the ejection velocity as 
a result of a Deep Impact-like impact. The same initial density (i.e., the same 
mass) was used for the porous and non porous cases. 

Figure 10: Impact in a basalt target with an initial distention of ao = 3.0 
(black). Full compaction (a = 1, white) is only reached using 1 model (top). 
Using model 2 (middle) or 3 (bottom), there are no fully compacted particles. 
For the definition of the models, see text. 

Figure 11: Distention as a function of distance (negative z-direction) obtained 
by model 1/2/3. 

Figure 12: Cumulative volume of ejected material for the models 1/2/3. For 
comparison, the total volume which was compacted is also shown for the three 
cases (horizontal lines). 

Figure 13: Cumulative volume of ejected material for model 1 using = 
2x10^ dyn/cm^ and Pg = 1x10^ dyn/cm^. The horizontal lines indicate the 
total volume which was compacted in each case. 
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Figure 1: 
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